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Abstract 

NASA and its contractors are working on 
structural concepts for absorbing impact energy of 
aerospace vehicles. Recently, concepts in the form of 
multi-cell honeycomb-like structures designed to 
crush under load have been investigated for both 
space and aeronautics applications. Efforts to 
understand these concepts are progressing from tests 
of individual cells to tests of systems with hundreds 
of cells. Because of fabrication irregularities, 
geometry irregularities, and material properties 
uncertainties, the problem of reconciling analytical 
models, in particular LS-DYNA models, with 
experimental data is a challenge. A first look at the 
correlation results between single cell load/deflection 
data with LS-DYNA predictions showed problems 
which prompted additional work in this area. This 
paper describes a computational approach that uses 
analysis of variance, deterministic sampling 
techniques, response surface modeling, and genetic 
optimization to reconcile test with analysis results. 
Analysis of variance provides a screening technique 
for selection of critical parameters used when 
reconciling test with analysis. In this study, complete 
ignorance of the parameter distribution is assumed 
and, therefore, the value of any parameter within the 
range that is computed using the optimization 
procedure is considered to be equally likely. Mean 
values from tests are matched against LS-DYNA 
solutions by minimizing the square error using a 
genetic optimization. The paper presents the 
computational methodology along with results 
obtained using this approach. 

Introduction 

The problem of updating parameters in finite 
element models using experimental data is one that 
has been studied for years. Although computer 
capability has increased dramatically and the size of 
models is perhaps at record highs, the task of 


reconciling a model with test has not changed. Because 
the approaches used are highly diverse, several 
engineering societies have produced validation and 
verification guides Ref. [1-2] to aid users through the 
process. In spite of these, each problem presents unique 
challenges and even the application of existing 
approaches often require development of new metrics. 
From the computational viewpoint, many advances have 
been made to address different aspects of the problem; 
for example, response surface techniques to approximate 
solutions for computationally intensive problems, 
deterministic sampling approaches to improve coverage 
of parameter variations, efficient optimization algorithms 
to handle the most difficult problems, design of 
experiments techniques to plan and sample the design 
space, and of course the modeling tools to handle the 
most complex problems. This work uses LS-DYNA, 
Ref. [3], as the primary modeling tool and even though 
LS-DYNA offers LS-OPT [4] for design optimization, 
the work here uses scripts and tools developed within 
MATLAB, (Copyrights 1984-2007, The Mathworks 
Inc). 

The example problem selected to demonstrate the 
update approach is that of a three point bend test of a 
single energy absorbing cell, loaded specifically to 
understand its behavior in shear. Of particular interest is 
the buckling or failure load of the cell since it relates 
directly to the efficiency of the honeycomb concept. 
Although this is a static loading problem, in this paper 
LS-DYNA [3] explicit solution approach is used with 
loads incrementally applied as a function of time until 
the full load is reached or the model fails. 
Load/displacement test results for the single cell model 
have been collected and the task is to recover a set of 
parameters to reconcile the test and the LS-DYNA 
model. 

To update parameters using experimental data, there 
are several key ingredients in the mathematical 
formulation that must be combined to make the approach 
computationally efficient: the use of response surface 
models specifically, the Moving Least Squares (MLS) 
[5] and Kriging’s method, [6-7], the use of a Halton-leap 


deterministic sampling technique [8-9], the use of 
analysis of variance for parameter screening [10], and 
the use of a genetic algorithm for optimization [11]. 
Admittedly, all these are well established approaches 
documented in the open literature, but the 
combination of the various elements for parameter 
update is not so commonly used, and that is the main 
thrust of this work. Aside from the algorithms, there 
are other important aspects of any model update 
effort that must be kept in mind when choosing a 
solution strategy. For example, the possibility of 
rejecting a model because it is inadequate should 
always be an acceptable outcome from any solution 
strategy. Also important is the reliance on 
engineering judgment when prescribing parameter 
variations to ensure feasible parameter sets. All these 
elements have been incorporated in the solution 
strategy discussed in the paper. 

Description of the LS-DYNA Model 

The testbed and LS-DYNA model used in this 
work were presented in Ref. [12] and so only a brief 
description is presented here for the reader to get 
familiar with the model. Figure 1 shows the model 
consisting of three major parts: the “Bakelite” 
stanchions, the hexagonal-shaped deformable 
composite cell-wall, and the flanges. The cell 
dimensions are cell-wall width (W) = 0.75”, cell-wall 
thickness (t) = 0.010” and overall length equal to 
3.35”. The model was constructed using an element 
edge length of 0.05”. The Kevlar cell-wall is 
modeled using *MAT_58 or 

*MAT_LAMINATED_COMPOSITE_FABRIC in 
the LS-DYNA model. Data from tensile tests of 
0°/90° coupons were used to determine approximate 
values of modulus and failure strain in both the 
longitudinal and the transverse directions. Likewise, 
data obtained from tensile tests of ±45° coupons were 
used to estimate the shear stiffness and strength. The 
flanges consisted of two layers of Kevlar fabric, with 
a total effective thickness of 0.02”. The thickness of 
the adhesive layer in the flanges was neglected. The 
specific nominal material properties used in the 
model to represent the Kevlar 129 fabric are shown in 
Table 1. The Bakelite stanchion was represented 
using a linear elastic material with a density of 
1.356E-4 lb-s 2 /in 4 , an elastic modulus of 1.09E+6 psi, 
and a Poisson’s ratio of 0.25. The complete model 
consisted of 27,582 nodes; 16,840 hexagonal solid 
elements; and 8,040 Belytschko-Tsay shell elements. 

The test specimen was loaded in three-point 
bending by applying a compressive load at the center 
of the hexagonal cell using a 0.25” wide bar, as 
shown in Fig. 1. In the test, loading was introduced 
with displacement control with a load rate of 2.0 ipm. 
The reaction points (0.25” from the bottom edge) 


were assumed to be simply supported. To represent the 
simply supported boundary conditions, a 
*BOUNDARY_SPC_SET was defined in the model. For 
this set, all of the nodes within the 0.25” by 0.75” area at 
the bottom of both ends of the specimen were 
constrained initially and later changed to a line 
constraint. 

Model Reconciliation Metric 

To reconcile test results with analysis, a metric must 
be established to quantify closeness of the LS-DYNA 
predictions with test. Because load versus deflection 
data is available, for a test displacement vector y (t) and 

the LS-DYNA response for a given parameter vector 
vTs y(t,v ) , the squared error difference is, 

^(K) = Z(T(^Ar,rt,.)-j,(Mr)) 2 (i) 

k-0 

where N is the number of data points and AT is the 
elapsed time between samples. Note that this metric is 
written as a sum of error differences over time to 
highlight the fact that the explicit solution in LS-DYNA 
is used. On the other hand, since the load is increased at 
every time step, load values are proportional to time; that 
is, if the maximum load during a test is P mlx and the 
total simulation time (arbitrarily selected) is T f = NAT , 
the load increment is P / N ■ This ensures that the full 

max 

load is reached at the final time unless failure occurs. 

With the metric defined in Eq. (1), the LS-DYNA 
parameters can now be adjusted to minimize the error 
difference over all the possible parameter sets and all 
experimental measurements. In this paper, a 
deterministic solution is sought, but if the parameter 
variations are prescribed in terms of distribution 
functions, it is straightforward to find the most likely 
parameter set along with its probability. 

Computational Framework 

To conduct a study like this it is preferable to 
automate the generation of LS-DYNA solutions and 
parameter changes. It is also important to develop 
computational tools that can take advantage of models as 
developed by engineers. Because the impact dynamics 
community uses LS-DYNA routinely, the most flexible 
approach is to manipulate LS-DYNA input file structure 
directly. Figure 2 shows a data flow diagram 
implemented using MATLAB Script files. These script 
files modify the LS-DYNA keyword input file 
automatically to update parameter values using a priori 
knowledge of the parameter distributions, execute LS- 
DYNA, and read LS-DYNA output files. By storing the 
results within the MATLAB environment, all the 
MATLAB toolboxes are available for use. 



In order to make this approach viable for 
computationally intensive LS-DYNA models, it is 
proposed (as depicted in the center of figure 2), that 
input-output mapping of the parameter values to LS- 
DYNA response outputs be captured using an 
adaptive response surface technique. For this task 
two critical elements are required; 1) an efficient 
response surface technique, and 2) an efficient multi- 
dimensional sampling technique. Comments on the 
selection of both approaches are provided next. 

Moving Least Squares (MLS) Response Surface 
Formulation 

A response surface model is a mathematical 
representation of input variables (variables that the 
user controls) and output variables (dependent 
variables). Many papers have been published on 
response surface techniques and one of the 
approaches selected was developed by 
Krishnamurthy as described in Ref. [5]. In this 
formulation the input/output relationship is given in 
parametric form as 

u = p t a-'bu 

A = 'Zw l (y)piv l )pf(v l ) ^ 

5 =Z w X v )^ v i) 

1=1 

pT = bO,) p(y 2 ) ••• /-On)] 
where U e M lx?iv is a vector of predictions, 
U g R 1 ^ is a vector of responses (in this case LS- 
DYNA solutions and stacked row-wise), yds the ; th 

parameter vector from a sample population whereas 
V is a variable representing the parameters, N is the 
population size (also number of time steps), w (v) is a 

user-defined function that weights the proximity of 
other parameter vectors on the response surface, q is 
the number of outputs (sensors), and p(v) is a set of 
basis functions. Krishnamurthy provided several 
weighting functions to handle problems with 
different continuity requirements given as a function 
of the proximity radius, where the radius was defined 
as p = II v - vll // and / is a user defined distance. In 

our implementation of MLS, the proximity radius is 
computed directly from data using a quadratic search 
to minimize the error between the data and the 
response surface prediction. Also to avoid having a 
catalog of weighting functions for problems with 
different continuity requirements, the 
sinc(/?) = sin(/?) / p function is used instead. 

To assess accuracy of the response surface 
technique from a multidimensional viewpoint, the dot 
product of the response surface prediction vector and 


the LS-DYNA solution vector is used to estimate the 
error in vector angle and magnitude. For instance, 
consider Y e R** 1 to be a vector of response surface 
predictions for q output nodes and Y e to be a 
vector of LS-DYNA predictions, the angle error is 

E g =1 -cos(0) =1-^4- ( 3 ) 

|t||t 


and the magnitude error is 


E m 



(4) 


These two metrics are used later to compare techniques. 


Kriging’s Response Surface 

This alternate response surface formulation has been 
thoroughly studied and documented in numerous papers; 
two excellent references are [6] and [7]. In this approach 
estimates of the response are computed using 

y{v) = P + >\v) T R-\Y-pi ) ( 5 ) 

with 

P = (f R-'lp'f R-'Y 

m 

R( v i’ v j) = exp[-0^T (vf - v*) 2 ] 

k=l 

for i=l,2,...,N and 

r T (y) = [R{v,v x ) R(v,v 2 ) ... R(v,vJ\ 

Y T =[)\ T 2 ••• Tvl 

where y and [j are scalars, R e R WxJV is a correlation 
matrix, r(v) e R Wxl is a correlation vector, Y g M :VxI is a 
vector of responses, 1 e R iVxl is a vector whose 
elements are all ones, 9 is a scalar parameter yet to be 
computed, and m is the number of design variables. 
Although Eq. (5) is written for a system with a scalar 
output, the formulation can be applied sequentially for 
cases with multiple outputs. A subtle but critical aspect 
of this formulation is that the unknown parameter 9 , 
used to define the correlation matrix R, must be 
computed using maximum likelihood estimation. That is 
9 is the value that maximizes the scalar metric 

Z(6») = -i(Aln(<r 2 ) + ln|R|), 6>g[0,oo] (6) 

For an estimated standard deviation 

G-=fY-p\) T R-\Y-p\) 

Solving this scalar optimization problem produces a 
response surface fitting error that is distributed according 
to a Gaussian distribution with zero error at the known 
solutions. In contrast to the MLS approach, computing 
the parameter 9 requires a scalar optimization for each 
output vector. Hence, if the problem at hand contains q 
outputs sampled N times to generate a Kriging’s surface 
one needs qN optimizations of the metric in Eq. (6). 



Although computationally intensive, the solutions 
can be very accurate as will be demonstrated later in 
the paper. 

Deterministic Sampling of the input Parameters 

A critical step when creating response surface 
models is in the sampling of the parameter domain. 
That is, having selected a set of parameters as our 
inputs (control variables) to the response surface 
algorithm, sampling of parameter values over their 
prescribed domain is critical. For this purpose the 
Halton-leap low discrepancy deterministic sampling 
approach described in Ref. [8] and studied 
extensively in Ref. [9] has been selected. The 
selection is based not only on the improved 
convergence of statistical parameters that this 
approach provides (over strictly random sampling) 
but also in that it allows setting of the problem 
sequentially. For example, if the number of LS- 
DYNA runs required is unknown a priori, one may 
need to add solutions later in the process. Using 
Flalton-leap, this procedure is simple to do without 
the risk of generating repeated solutions or sacrificing 
coverage, which is not the case with many 
conventional sampling techniques. 

Discussion of Results 

The model as developed and reported in Ref. 
[12] is initially used to assess the effects of parameter 
uncertainties in the predicted displacements. At first, 
engineering judgment is used to select a set of 
parameters that is deemed uncertain. In addition, 
parameters that control the execution of LS-DYNA 
are also included. Of course, these parameters can 
appear anywhere in the LS-DYNA input file but for 
our purposes, only those in the MAT 58 are 
considered. The top seven parameters in Table 1 are 
chosen and parameter values are generated using 
Flalton-leap deterministic sampling. Initially, 40 LS- 
DYNA runs are executed to create a 2 nd order 
response surface (admittedly only 36 solutions are 
needed) and to conduct analysis of variance. 

Nominal Test and LS-DYNA comparisons 

Figure 3a shows the load versus displacement 
variations as parameters are varied over their pre- 
defined ranges. The response envelope; test results 
upper bound, lower bound, and mean values from 
five experiments are plotted versus predicted upper, 
lower, and mean values from LS-DYNA; for a cell 
model with bottom edge area constraint. Bounds for 
the analysis data are computed from results after 
arbitrarily vaiying all parameters over their domain. 
After examining Fig. 3a it is apparent that the LS- 
DYNA response envelop, for arbitrary parameter 
variations, does not bound the measured values. This 


prompted another look at the model that resulted in 
changes to the end constraint from an area to a line 
constraint. Because the initial stiffness is also over 
predicted, parameter 8 is added to allow for variations in 
the longitudinal and transverse modulus. Figure 3b 
shows the LS-DYNA predictions after these changes 
were made showing that now a large portion of the test 
data envelop is bounded by the LS-DYNA predictions. 

Analysis of variance 

In studies like this, the first step is always parameter 
importance screening, which is often conducted using 
analysis of variance. For this, the conventional approach 
is to use the design of experiment (DOE) methodology 
[10] to prescribe variations of the controlled variables 
(model parameters) and to compute the variance effects 
directly from the observations. Here, a deterministic 
sampling approach is used instead to prescribe parameter 
variations, and variable screening is conducted directly 
from the response surface by varying each parameter 
independently. Although this approach sacrifices some of 
the advantages of DOE, coverage of the parameter space 
is thought to be better for problems where the goal is to 
create accurate response surface models as opposed to 
minimize the number of experiments. Certainly, if the 
surface is not accurate, the variance estimate is in error, 
but for parameter screening purposes this approach is 
considered adequate. Results from the analysis of 
variance are shown in Figure 4 for the parameters in 
Table 1. The abscissa of Fig. 4a indicates parameter 
numbers ranging from 1 to 8 and the ordinate shows five 
outputs corresponding to nodes 529-533 in the LS- 
DYNA model (located along the area where the loads are 
applied). Because the variance analysis is conducted at 
each load increment (or time increment), Fig. 4a shows 
the case for 12.4 lbs first. When all the parameters are 
varied simultaneously the displacement at each node will 
vary and the maximum standard deviation a over all 

the parameter variations and all outputs is computed and 
used to normalize the data. At a particular output, the 
ratio of its displacement standard deviation a to cr max is 

a measure of the how much displacement occurs at that 
output node. Figure 4a shows cr/cr with blue squares 

sized according to the numerical value of the ratio (hill 
square cr/<T max = 1 )• Because the load is symmetric and 

all nodes deform practically the same amount. Fig. 4a 
shows full size squares throughout. For cases where the 
load is not symmetric, one would see variations in this 
ratio as a function of the output location. This procedure 
helps to highlight those nodes with maximum 
displacements. 

Now for parameter screening, it is necessary to 
determine how much a particular parameter variation 
contributes to the total displacement variance. By 



varying one parameter at a time, the proportion of the 
total variance contributed by variations of a particular 
parameter is computed and superimposed on Fig. 4a 
using squares with black lines. Thus, Fig. 4a shows 
that at a load of 12.4 lbs parameter 3 (shear modulus) 
and 8 (longitudinal and transverse modulus) are the 
two largest contributors to the displacement variance. 
In addition to the variance, one should also consider 
the effects of parameter variations on the mean 
response. For this case, the mean value 
corresponding to each output in printed to the right of 
Fig. 4a. Also the center of the black squares 
coincides with the center of the blue squares only 
when both mean values (that is the output mean and 
the mean from parameter changes), are the same. 

In contrast, Fig. 4b shows the same variance 
information but for a load case of 2741bs. Note that 
the variance contribution from parameters 3 and 8 is 
now insignificant as compared to the contribution 
from parameters 4 (SLIMS) and 5 (ERODS). 
Moreover, parameters 1, 6, and 7 also have 
noticeable contributions, but from Figs. 4a and 4b, 
parameter 2 (strain limit) had a very small impact on 
the displacement variations and thus could have been 
removed (however, it was not removed). An 
important conclusion drawn from the variance 
analysis is that parameter sensitivity (in terms of 
variance) is a function of load (time) and must be 
computed as such or risk removing critical variables 
from the problem. 

Evaluation of Response Surface Techniques 

With LS-DYNA solutions at hand and response 
surfaces computed using both MLS and Kriging, the 
phase E g and magnitude error metrics E m defined in 

Eqs. (3-4) are easily computed at each known LS- 
DYNA solution. Figure 5a shows E using MLS and 

E g as a function of the solution number, which 

ranged from 1 to 87. Likewise, on Fig. 5b are the 
corresponding results using Kriging’s approach. 
Note that the Kriging solution is orders of magnitude 
more accurate than MLS when using this metric 
because it forces exact matching of LS-DYNA 
solutions. This accuracy is only possible after 
optimizing for a load/time varying solution for 9 . 

Model Reconciliation and Optimization 

After a satisfactory response surface model is 
initially created, the next step is to use it to solve for 
a parameter set that minimizes the performance 
metric in Eq. (1). This optimization process is 
iterative in nature. Figure 6 shows intermediate 
results for cycle 2 of the optimization process. 
Specifically, Fig. 6a shows the MLS predictions, the 
test results envelope, and LS-DYNA predictions 


when using cycle 2 optimized parameters, whereas the 
plot in Fig. 6b shows Kriging’s results. In the context 
of this work, a cycle (or iteration) starts by creating a 
response surface model from an initial set of LS-DYNA 
solutions, as depicted in Fig. 7, followed by an 
estimation of an optimum parameter set. To verify this 
optimum solution, LS-DYNA is executed using the 
optimized parameter set and this solution is compared to 
our response surface predictions. This comparison is 
shown in Fig. 6 after going through the cycle twice. If 
the comparison is satisfactory the optimization cycle can 
be stopped. Otherwise, this suboptimal LS-DYNA 
solution is added to the LS-DYNA solution set and an 
updated response surface model is created. At this point, 
the cycle starts again until convergence is achieved. 
Naturally, solutions after each cycle with MLS and 
Kriging are different because the response surface 
models are different. For example, Table 2 shows the 
optimized parameter set obtained using MLS and 
Kriging’s with 155 LS-DYNA solutions. Figure 8 shows 
the results using Kriging’s method after completing 3 
cycles. The plot in figure 8a shows mean values for 
Kriging, LS-DYNA, and the mean of the test results. The 
plot shown in figure 8b includes the same data but the 
test results are shown in terms of bounds. Because 
Kriging and LS-DYNA predictions are close to the mean 
test results, the optimization cycle was stopped. 
Conversely, the MLS optimization did not converge after 
4 cycles and so results are not shown. 

Concluding Remarks 

A methodology to update LS-DYNA generated 
models has been successfully demonstrated that uses 
response surface techniques to reduce computational 
burden, analysis of variance for parameter screening, a 
genetic optimization algorithm for parameter updates, 
and deterministic sampling for efficient sampling of the 
parameter space. To illustrate the approach, an LS- 
DYNA model of an energy absorbing composite cell is 
used and 8 parameters were updated based on 
experimental data. In this example, load/displacement 
data from five experiments are used to recover a 
parameter set to reconcile test with analysis. Initial 
analysis of variance was instrumental in understanding 
deficiencies in the modeling approach and pointed to 
improvements that resulted in better correlation even 
without optimization. Results show that the prediction 
error is significantly reduced and in the end, the LS- 
DYNA model accurately predicted the mean response 
measured experimentally. 

Because complete ignorance of the parameter 
uncertainties is assumed, it is not possible to assess the 
probability that this optimized parameter set is correct. 
Nonetheless, the computational framework can 
incorporate distribution functions that would allow users 
to estimate this probability as well. 



Finally, of the two different response surface 
techniques, Moving Least Squares (MLS) and 
Kriging, results using Kriging’s approach proved to 
be more accurate than MLS but the computational 
burden is significantly higher. 
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Table 1. MAT 58 Material Properties and distribution values 


Parameter 

No. 

Material Property 

Upper 

Bound 

Lower 

Bound 

Nominal 

1 

Stress limit of nonlinear portion 
of shear curve psi 

8000 

2000 

4156.3 

2 

Strain limit of nonlinear portion 
of shear curve in/in 

0.0374 

0.00935 

0.022509 

3 

Shear modulus, psi 

614,962.0 

1.5374e+5 

4.563e+5 

4 

SLIMS 

0.75 

0.1 

0.30087 

5 

ERODS 

0.8 

0.05 

0.28554 

6 

Strain at shear strength, in/in 

0.1 

0.01 

0.053669 

7 

Shear strength, psi 

14000.0 

3750 

8077 

8 

Young’s modulus, longitudinal 
and transverse, psi 

1.6E+6 

1.3e+006 

1.3183e+6 


Table 2. Optimized parameter values using MLS and Kriging 


Parameter 

No. 

Material Property 

Nominal 

Optimum 

MLS 

(155) 

Optimum 

KRG 

(155) 

1 

Stress limit of nonlinear 
portion of shear curve psi 

4156.3 

8000 

7020.8 

2 

Strain limit of nonlinear 
portion of shear curve in/in 

0.022509 

0.034969 

0.00935 

3 

Shear modulus, psi 

4.563e+5 

1.5374e+5 

1.5374e+5 

4 

SLIMS 

0.30087 

0.64059 

0.1 

5 

ERODS 

0.28554 

0.05 

0.73279 

6 

Strain at shear strength, in/in 

0.053669 

0.041606 

0.094463 

7 

Shear strength, psi 

8077 

3750 

12383 

8 

Young’s modulus 
longitudinal direction and 
transverse, psi 

1.3183e+6 

1.3e+006 

1.3e+006 








Deformable composite 



Fig. 1 Picture of cell testbed and LS-DYNA model 



Fig. 2 Computational framework for update problem 









a) Area constraint b) Line constraint 

Fig. 3. Comparison of original LS-DYNA model (3a) with area constraint versus 
model with line constraint (3b) 
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b) Load 274 lbs. Parameter No. 
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Fig. 4 Analysis of variance at two different loads for outputs 1-5 corresponding 
to displacements in the y direction for nodes 529-533 
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a) MLS 


b) Kriging 


Fig. 5. Comparison of displacement prediction error with MLS and Kriging. 



Fig. 6 Comparison of test results against MLS and Kriging’s response surface procedure with the 
optimized LS-DYNA solution (cycle 2) 
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Fig. 7 Optimization design cycle 



Fig. 8 Comparison of test mean and bounds against Kriging’s response surface model and 
LS-DYNA using the optimized parameters (cycle 3) 


